Aplicações de Modelos Lineares Generalizados
Elton Dantas de Oliveira Mesquita
Pacotes
library(tidyverse)
library(psych)
library(plotly)
library(GGally)
library(car)I. Modelo Linear Normal
Para a aplicação a seguir, consideraremos o modelo de regressão normal linear
\[ y_i = \beta_1 + \beta_2x_{2i} + ... + \beta_px_{pi} + \epsilon_i, \quad i = 1,...,n \]
em que os erros \(e_i\) são variáveis aleatórias independentes, normalmente distribuídas, de média zero e variância \(\sigma^2\) constante.
Previsão do Preço de Venda de Imóveis
Dados
Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 112). Dados disponíveis em: https://www.ime.usp.br/~giapaula/textoregressao.htm
Neste exemplo, vamos modelar o preço de venda de imóveis a partir de dados relativos a uma amostra de 27 imóveis. As variáveis do conjunto de dados são:
imposto: imposto do imóvel (em US$ 100);
areaT: área do terreno (em 1.000 pés quadrados);
areaC: área construída (em 1.000 pés quadrados);
idade: idade da residência (em anos);
preco: preço de venda do imóvel (em US$ 1.000).
Sendo assim, o nosso objetivo é encontrar o melhor ajuste linear que explique e quantifique a variável preço de venda a partir das demais.
# Obtendo os dados
dados = read.table("dados/imoveis.dat")
colnames(dados) = c("imposto", "areaT", "areaC", "idade", "preco")
attach(dados)
# Algumas observações dos dados
head(dados)## imposto areaT areaC idade preco
## 1 4.9176 3.472 0.998 42 25.9
## 2 5.0208 3.531 1.500 62 29.5
## 3 4.5429 2.275 1.175 40 27.9
## 4 4.5573 4.050 1.232 54 25.9
## 5 5.0597 4.455 1.121 42 29.9
## 6 3.8910 4.455 0.988 56 29.9# Medidas descritivas
describe(dados)## vars n mean sd median trimmed mad min max range skew
## imposto 1 27 7.24 2.88 6.09 6.84 2.15 3.89 15.42 11.53 1.40
## areaT 2 27 6.35 2.40 5.85 6.22 2.07 2.28 12.80 10.53 0.68
## areaC 3 27 1.51 0.56 1.49 1.41 0.39 0.98 3.42 2.44 2.02
## idade 4 27 36.48 14.05 40.00 36.96 14.83 3.00 62.00 59.00 -0.34
## preco 5 27 38.50 14.31 36.90 35.65 8.90 25.90 84.90 59.00 2.25
## kurtosis se
## imposto 1.34 0.55
## areaT 0.01 0.46
## areaC 4.08 0.11
## idade -0.54 2.70
## preco 4.63 2.75Análise Descritiva
Box-plots
plot_ly(type = "box") %>%
add_trace(y = imposto, name = "imposto") %>%
add_trace(y = areaT, name = "área do terreno") %>%
add_trace(y = areaC, name = "área construída") %>%
add_trace(y = idade, name = "idade do imóvel") %>%
add_trace(y = preco, name = "preço de venda")Em um breve resumo descritivo dos dados, observamos que as altas variâncias das variáveis idade e preço destacam-se das demais. Com uma idade mediana de 40 anos, temos um imóvel de apenas 3 anos de idade e, com um preço mediano de US$ 36.900, temos dois imóveis bem mais caros nos valores de US$ 82.900 e US$ 84.900.
Os gráficos de box-plot nos mostram que os dados possuem outliers e que as variáveis possuem distribuições assimétricas. Focando na variável alvo, preço de venda, os pontos destoantes são justamente os dois imóveis mais caros.
Densidades, Dispersões e Correlações
g = ggpairs(dados, aes(color = I("slategray"), fill = I("slategray")),
lower = list(continuous = wrap("smooth",col="black")),
diag=list(continuous = wrap("densityDiag",alpha=0.5,size=1)))
ggplotly(g)No gráfico acima, as curvas representadas na diagonal principal deixam mais evidente a constatação de assimetria nas distribuições observada nos box-plots. A variável preço possui correlações positivas fortes com as variáveis imposto, área do terreno e área construída, mas uma correlação negativa fraca com a variável idade. Em outras palavras, quanto maiores os valores de imposto, área construída e área do terreno, maior o preço de venda do imóvel.
Podemos observar também que há correlações relevantes entre as variáveis explicativas imposto, área construída e área do terreno. Isso nos dá indícios de multicolinearidade e, assim, uma possível redundância de informação passada por essas variáveis. Então, talvez um modelo completo com todas as variáveis não seja o mais adequado.
Modelo Completo
Nesse primeiro ajuste, consideraremos o modelo completo com todas as variáveis disponíveis.
\[ preço_i = \beta_0 + \beta_1imposto_i + \beta_2areaC_i + \beta_3areaT_i + \beta_4idade_i \]
model1 = lm(preco~.,dados)
summary(model1)##
## Call:
## lm(formula = preco ~ ., data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9685 -2.7152 0.2663 2.1374 6.9872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.43587 4.09233 0.595 0.55776
## imposto 2.07823 0.55296 3.758 0.00109 **
## areaT 0.23238 0.50658 0.459 0.65094
## areaC 13.97381 2.90663 4.808 8.4e-05 ***
## idade -0.04376 0.06628 -0.660 0.51594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.077 on 22 degrees of freedom
## Multiple R-squared: 0.9313, Adjusted R-squared: 0.9188
## F-statistic: 74.56 on 4 and 22 DF, p-value: 1.809e-12Verificamos que, para o modelo completo, as variáveis imposto e área construída são as únicas significativas e obtivemos um \(R^2\) ajustado de 0.92, muito alto, indicando um possível bom ajuste aos dados amostrais. Porém, o \(R^2\) por si só não nos garante que o modelo seja o mais adequado.
Seleção de variáveis
Aplicaremos o critério de informação de Akaike (AIC) ao modelo completo, para selecionar as variáveis que deverão permanecer.
MASS::stepAIC(model1)## Start: AIC=80.36
## preco ~ imposto + areaT + areaC + idade
##
## Df Sum of Sq RSS AIC
## - areaT 1 3.50 369.14 78.614
## - idade 1 7.25 372.88 78.887
## <none> 365.64 80.357
## - imposto 1 234.76 600.40 91.748
## - areaC 1 384.13 749.77 97.746
##
## Step: AIC=78.61
## preco ~ imposto + areaC + idade
##
## Df Sum of Sq RSS AIC
## - idade 1 11.51 380.64 77.443
## <none> 369.14 78.614
## - imposto 1 246.10 615.24 90.407
## - areaC 1 489.34 858.47 99.402
##
## Step: AIC=77.44
## preco ~ imposto + areaC
##
## Df Sum of Sq RSS AIC
## <none> 380.64 77.443
## - imposto 1 350.01 730.65 93.049
## - areaC 1 483.16 863.80 97.569##
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
##
## Coefficients:
## (Intercept) imposto areaC
## 0.7903 2.2971 13.9333O AIC nos retornou as variáveis imposto e área construída como aquelas que deverão ser mantidas no modelo, o qual chamaremos de modelo parcimonioso.
Modelo Parcimonioso
model2 = lm(preco~imposto+areaC,dados)
summary(model2)##
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0177 -3.3169 0.0958 2.4382 7.0949
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7903 2.2794 0.347 0.732
## imposto 2.2971 0.4890 4.698 8.96e-05 ***
## areaC 13.9333 2.5244 5.519 1.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.982 on 24 degrees of freedom
## Multiple R-squared: 0.9285, Adjusted R-squared: 0.9225
## F-statistic: 155.8 on 2 and 24 DF, p-value: 1.791e-14E, assim, obtemos um modelo em que as variáveis imposto e área construída são altamente significativas.
Porém, na análise descritiva, observamos que as variáveis área construída e área do terreno possuem uma correlação considerável. Então, podemos ainda testar um segundo modelo parcimonioso considerando a variável área do terreno ao invés da área construída.
model3 = lm(preco~imposto+areaT,dados)
summary(model3)##
## Call:
## lm(formula = preco ~ imposto + areaT, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.621 -3.070 0.446 3.167 10.610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1614 3.3019 0.957 0.3479
## imposto 3.9126 0.5293 7.391 1.24e-07 ***
## areaT 1.1017 0.6347 1.736 0.0954 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.655 on 24 degrees of freedom
## Multiple R-squared: 0.8558, Adjusted R-squared: 0.8438
## F-statistic: 71.22 on 2 and 24 DF, p-value: 8.081e-11E, assim, obtemos um modelo em que a variável área do terreno passa de não significativa para pouco significativa ao nível de 10% de significância. Já a variável imposto continua bastante significativa ao nível de 1% de significância.
Diagnóstico
Agora, faremos a análise de diagnóstico dos dois modelos parcimoniosos propostos, para verificar suas adequações.
Envelope
par(mfrow=c(1,2))
qqPlot(model2, pch = 16, main = "preco ~ imposto + areaC")## [1] 6 10qqPlot(model3, pch = 16, main = "preco ~ imposto + areaT")
## [1] 9 27Teste de Normalidade
\(H_0\): Os resíduos têm distribuição normal.
# Shapiro-Wilk
shapiro.test(model2$residuals)##
## Shapiro-Wilk normality test
##
## data: model2$residuals
## W = 0.96558, p-value = 0.4903shapiro.test(model3$residuals)##
## Shapiro-Wilk normality test
##
## data: model3$residuals
## W = 0.9793, p-value = 0.8462Ao nível de 5% de significância, não rejeitamos a hipótese de que os resíduos tenham distribuição normal. Nos gráficos de envelope para o ajuste normal, os pontos estão dentro das bandas, exceto pela observação 27 no modelo com área do terreno. Porém, a presença de leves ondulações pode ser indício de variância não constante. Vamos verificar se de fato isso ocorre no gráfico de resíduos por valores ajustados.
Resíduos x Valores ajustados
residualPlot(model2, type = "rstudent", id = TRUE,
ylab = "resíduo studentizado", xlab = "valor ajustado",
main = "preco ~ imposto + areaC", pch = 16, quadratic = FALSE)
residualPlot(model3, type = "rstudent", id = TRUE,
ylab = "resíduo studentizado", xlab = "valor ajustado",
main = "preco ~ imposto + areaT", pch = 16, quadratic = FALSE)
Embora não haja violação do pressuposto de normalidade dos resíduos, o gráfico de resíduos por valores ajustados nos dá indícios de heterocedasticidade, isto é, variância não constante, como havíamos suposto a partir do gráfico de envelope. Porém, esse comportamento pode estar sendo influenciado pelas observações destoantes vistas na etapa descritiva. Além disso, podemos perceber no gráfico de resíduos por valores ajustados a presença de pontos aberrantes, isto é, muito além do intervalo [-2,2]. No caso do primeiro modelo parcimonioso, o ponto aberrante é a observação 10, enquanto que no segundo modelo parcimonioso são as observações 9 e 27. Sendo assim, a seguir, vamos verificar se de fato esses pontos estão influenciando na qualidade dos ajustes.
Pontos Influentes
Partiremos agora para a investigação de possíveis pontos influentes, que podem estar atrapalhando a qualidade do ajuste.
Medidas de influência:
DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo;
DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores ajustados do modelo;
COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo;
HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados, é a métrica de alavancagem.
inf = influence.measures(model2)
summary(inf)## Potentially influential observations of
## lm(formula = preco ~ imposto + areaC, data = dados) :
##
## dfb.1_ dfb.imps dfb.areC dffit cov.r cook.d hat
## 9 -0.28 0.00 0.19 0.35 2.18_* 0.04 0.49_*
## 10 -1.25_* 0.29 0.59 1.62_* 0.87 0.73 0.32
## 27 -0.13 -2.04_* 1.86_* -2.12_* 2.04_* 1.39_* 0.61_*inf = influence.measures(model3)
summary(inf)## Potentially influential observations of
## lm(formula = preco ~ imposto + areaT, data = dados) :
##
## dfb.1_ dfb.imps dfb.areT dffit cov.r cook.d hat
## 9 -1.15_* 1.65_* -0.46 2.00_* 0.81 1.07_* 0.37_*
## 10 -1.21_* 0.47 0.68 1.54_* 1.02 0.69 0.35_*
## 27 0.09 -2.83_* 2.36_* -3.05_* 0.33_* 1.85_* 0.35_*Baseando-se nas medidas de influência acima, que verificam o efeito da remoção de cada uma das observações individualmente, temos que os possíveis pontos de influência são as observações 9 e 27, para ambos os modelos parcimoniosos, e a observação 10, para o segundo. Percebemos também que a observação 27 é a mais crítica, por impactar nas estimativas dos dois parâmetros, nos valores ajustados e na variância do modelo. Também podemos constatar isso graficamente, como se segue a baixo.
# DISTANCIA DE COOK
plot(model2, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
sub.caption = "preco ~ imposto + areaC")
plot(model3, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
sub.caption = "preco ~ imposto + areaT")
# ALAVANCAGEM
cut = 2*3/27 # 2*p/n
plot(hatvalues(model2), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem",
sub = "preco ~ imposto + areaC")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model2) * as.integer(hatvalues(model2) > cut),
cex=0.8, p=1, offset=0.3)
cut = 2*3/27 # 2*p/n
plot(hatvalues(model3), pch = 16,
xlab = "observação", ylab = "medida h", main = "Alavancagem",
sub = "preco ~ imposto + areaT")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model3) * as.integer(hatvalues(model3) > cut),
cex=0.8, p=1, offset=0.3)
As observações 9, 10 e 27 são referentes aos dados:
dados[c(9,10,27),c("imposto","areaC","areaT","preco")]## imposto areaC areaT preco
## 9 15.4202 3.42 9.8 84.9
## 10 14.4598 3.00 12.8 82.9
## 27 12.0000 1.20 5.0 41.0Os imóveis 9 e 10 possuem os valores de imposto mais elevados e as maiores áreas de terreno e de construção. De acordo com a relação linear identificada na etapa descritiva e a significância dessas variáveis para o modelo preditivo, podemos dizer que essas características justificam um alto preço de venda. Porém, os valores atribuídos aos imóveis 9 e 10 são muito mais elevados que os esperados pelo modelo proposto a imóveis com os mesmos atributos. Já a observação 27 causa estranheza por ter um valor de imposto tão alto e próximo aos dos imóveis 9 e 10, mesmo tendo menos da metade de suas áreas de terreno e de construção e custando pouco menos da metade que o imóvel 10.
Impacto das remoções das observações
# Compara os coeficientes e seus respectivos p-valores de um modelo
# após a remoção de uma observação
impacto = function(m1, m2){
imp = (coef(m1) - coef(m2))/coef(m1) * 100
impacto = data.frame(coef_com_i = coef(m1),
coef_sem_i = coef(m2),
impacto = imp,
pval_com_i = summary(m1)$coefficients[,4],
pval_sem_i = summary(m2)$coefficients[,4])
return(impacto)
}Observação i=9:
model2_sem9 = lm(preco~imposto+areaC,dados,subset=-9) impacto(model2,model2_sem9)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i ## (Intercept) 0.7903027 1.433603 -81.39918902 7.318305e-01 0.6304908219 ## imposto 2.2970580 2.297773 -0.03113799 8.955810e-05 0.0001221436 ## areaC 13.9332510 13.454944 3.43284823 1.122859e-05 0.0001144541model3_sem9 = lm(preco~imposto+areaT,dados,subset=-9) impacto(model3,model3_sem9)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i ## (Intercept) 3.161412 6.558508 -107.45502 3.478823e-01 5.383568e-02 ## imposto 3.912561 3.131158 19.97165 1.243498e-07 1.077092e-05 ## areaT 1.101699 1.360792 -23.51762 9.543188e-02 2.722222e-02A remoção da observação 9 não teve impactos relevantes nos coeficientes dos dois modelos parcimoniosos e também não causou mudança inferencial, uma vez que as variáveis imposto, área construída e área do terreno permaneceram igualmente significativas nos respectivos modelos;
Observação i=10:
model2_sem10 = lm(preco~imposto+areaC,dados,subset=-10) impacto(model2,model2_sem10)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i ## (Intercept) 0.7903027 3.398847 -330.069021 7.318305e-01 1.640664e-01 ## imposto 2.2970580 2.167849 5.624971 8.955810e-05 7.657851e-05 ## areaC 13.9332510 12.571459 9.773685 1.122859e-05 2.390378e-05model3_sem10 = lm(preco~imposto+areaT,dados,subset=-10) impacto(model3,model3_sem10)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i ## (Intercept) 3.161412 6.8915170 -117.988556 3.478823e-01 6.463495e-02 ## imposto 3.912561 3.6816714 5.901246 1.243498e-07 2.123188e-07 ## areaT 1.101699 0.6967438 36.757341 9.543188e-02 2.749152e-01A remoção da observação 10 não causou impactos relevantes nos coeficientes do primeiro modelo parcimonioso e também não ocasionou mudança inferencial, tendo as variáveis imposto e área construída mantido seus níveis de significância. Já no segundo modelo parcimonioso, a remoção da observação 10 ocasionou uma mudança inferencial, uma vez que a variável área de terreno perde sua significância;
Observação i=27:
model2_sem27 = lm(preco~imposto+areaC,dados,subset=-27) impacto(model2,model2_sem27)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i ## (Intercept) 0.7903027 1.068075 -35.14754 7.318305e-01 0.6320276693 ## imposto 2.2970580 3.255657 -41.73159 8.955810e-05 0.0001915344 ## areaC 13.9332510 9.412139 32.44836 1.122859e-05 0.0155595300model3_sem27 = lm(preco~imposto+areaT,dados,subset=-27) impacto(model3,model3_sem27)## coef_com_i coef_sem_i impacto pval_com_i pval_sem_i ## (Intercept) 3.161412 2.9423870 6.928084 3.478823e-01 2.603227e-01 ## imposto 3.912561 5.0694489 -29.568550 1.243498e-07 4.765596e-10 ## areaT 1.101699 -0.0528502 104.797154 9.543188e-02 9.260592e-01A remoção da observação 27 impactou em mudança inferencial nos dois modelos parcimoniosos. Enquanto que no primeiro, a significância da variável área construída é reduzida, no segundo, a variável área do terreno perde completamente sua significância.
Como o segundo modelo parcimonioso não nos trouxe nenhuma melhoria para o ajuste e as conclusões sobre os pontos de influência não favoreceram a permanência da variável área do terreno, optaremos pelo primeiro modelo parcimonioso com as variáveis imposto e área construída sugerido inicialmente pelo critério de seleção de Akaike.
A partir dos gráficos de resíduos por valores ajustados, de distância de Cook e de alavancagem, podemos fazer as seguintes classificações:
Ponto aberrante: observação 10;
Pontos de alavanca: observações 9, 10 e 27;
Ponto influente: observação 27.
Interpretação
Então, concluímos que
\[ preço_i = 0.79 + 2.30*imposto_i + 13.93*areaC_i \quad , \]
onde:
A cada US$ 100,00 acrescidos no valor de imposto, aumenta-se, em média, US$ 2.300,00 no preço de venda;
A cada 1.000 metros quadrados de área construída acrescidos, aumenta-se, em média, US$ 13.930,00 no preço de venda.